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Abstract 


This  report  presents  the  theoretical  prediction  of  the  low  frequency  un¬ 
steady  broadband  thrust  for  propellers/rotors  in  a  turbulent  flow,  such 
as  those  on  any  surface  ship.  Broadband  thrust  forces  are  predicted  by 
the  use  of  a  new  correlation  method  and  are  compared  with  existing  the¬ 
oretical  models  based  on  both  correlation  and  spectrum  approaches.  In 
contrast  to  the  available  theories,  the  new  theory  predicts  low  blade  rate 
humps  and  broadband  forces  simultaneously.  For  a  homogeneous  and 
isotropic  turbulent  inflow  and  for  a  given  advance  coefficient,  turbulence 
length  scale  and  turbulence  level,  the  present  theory  produces  an  unsteady 
broadband  thrust  level  which  is  proportional  in  the  high  frequency  region 
to  the  third  power  of  ship  speed.  At  low  frequencies,  predictions  of  pro¬ 
peller  unsteady  thrust  spectra  show  peaks  near  the  first  and  second  blade 
rate  frequencies  and  those  peaks  skew  to  the  higher  frequency  side  of  the 
blade  rate  frequencies.  The  physics  of  the  phenomena  are  discussed  and 
numerical  results  are  compared  with  limited  experimental  data. 

Administrative  Information 

This  work  was  performed  at  the  David  Taylor  Research  Center  (DTRC) 
in  Bethesda,  Maryland  20084-5000  with  work  unit  1-1904-301,  with  fund¬ 
ing  from  the  Office  of  Naval  Research,  element  0101224N. 

Introduction 

The  unsteady  forces  generated  by  a  propeller/rotor  consist  of  periodic 
force  components  and  a  broadband  force  component.  The  narrowband  pe¬ 
riodic  force  components  occur  at  multiples  of  the  blade  passage  frequency 
and  are  considered  to  result  from  the  unsteady  pressure  distribution  on  the 
propeller  blades.  These  unsteady  pressures  are  caused  by  the  non-uniform 
inflow  or  by  cavitation,  such  as  for  the  propellers  on  surface  ships.  Force 
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spectra  from  this  non-uniform  inflow  are  narrowband  at  the  blade  rates. 
The  high  frequency  broadband  noise  is  due  to  vortex  shedding  by  the 
propeller  blades.  Although  designers  are  careful  to  reduce  these  unsteady 
force  sources,  laboratory  experiments  indicate  that  an  unexpected  broad¬ 
band  force  occurs  at  first  and  second  blade  rate  frequencies.  This  low 
frequency  broadband  force  is  due  to  the  interaction  of  the  inflow  turbu¬ 
lence  with  the  propeller  blades.  The  inflow  turbulence  ingested  on  the 
propeller  blades  is  generated  in  the  hull  boundary  layer  and  in  the  bound¬ 
ary  layers  on  all  the  appendages  upstream  of  the  propeller,  superposed  on 
the  ambient  free-stream  turbulence.  Therefore,  the  inflow  turbulence  mav 
contain  a  broad  range  of  length  scales.  This  report  addresses  the  effects  of 
turbulence  on  broadband  forces  generated  by  propeller  blades  regardless 
of  the  source  of  the  turbulence. 

The  unsteady  force  caused  by  turbulence  ingestion  has  been  studied  by 
several  investigators  in  the  past.  Most  previous  studies  are  based  upon 
special  assumptions  and,  thus,  lack  generality.  Sevik  [1]  considered  the  for¬ 
ward  speed  to  be  much  greater  than  the  propeller  rotational  speed  and, 
therefore,  regarded  the  blades  as  rotationally  stationary  with  respect  to 
the  turbulence.  This  assumption  implied  that  the  blades  are  rotationally 
uncorrelated.  He  measured  the  unsteady  thrust  on  a  ten-bladed  rotor 
operating  in  a  water  tunnel  downstream  of  screen-generated  turbulence. 
The  experimental  measurements  show  humps  near  blade  rate  frequencies 
and  its  multiples  which  were  not  predicted  from  Sevik ’s  theory.  Thomp¬ 
son  [2]  extended  this  theory  and  performed  additional  experiments  for 
rotors  with  different  numbers  of  blades.  The  experimental  results  show 
that  the  magnitude  of  the  broadband  humps  increases  when  the  blade-to- 
blade  spacing  decreases.  He  included  in  his  analysis  the  blade  to  blade 
circumferential  correlations  in  a  zero  forward  speed  sense.  In  comparison 
to  the  experiments,  his  theoretical  results  produced  consistently  higher 
humps.  Chandrashekhara  [3]  examined  only  the  cases  for  which  the  cir¬ 
cumferential  correlation  length  was  larger  than  a  blade  spacing.  Mani 
[4],  Homicz  and  George  [5],  Amiet  [6],  and  Blake  [7]  considered  the  ro- 
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tor  forces  and  noises  radiated  by  turbulence  with  a  given  wave  spectrum. 
Mani  [4]  showed  that  the  radiated  sound  spectrum  from  a  rotor  operating 
in  a  turbulent  inflow  has  humps  centered  at  the  blade  rate  frequencies 
and  its  multiples.  The  width  of  those  humps  was  found  in  his  experi¬ 
ments  to  be  related  to  the  ratio  of  the  turbulence  length  scale  to  the  blade 
spacing.  Blake  [7]  considered  the  case  of  a  rotor  with  a  small  advance 
coefficient  and  a  turbulence  length  scale  considerably  smaller  than  the 
length  of  the  blade  span.  Simple  results  are  obtained  from  asymptotic 
expansion  of  small  and  large  correlation  lengths.  The  small  correlation 
result  approaches  the  result  of  Sevik’s  analysis,  and  the  large  correlation 
theory  predicts  the  narrow  bandwidth  humps  at  blade  rate  frequencies. 
Although  Blake’s  model  predicts  the  humps  at  the  blade  rate  frequency, 
the  hump  amplitudes  above  the  smooth  broadband  force  spectrum  curve 
do  not  exhibit  the  decay  phenomena  observed  in  experiments.  Recently, 
Martinez  [8]  extended  Sevik’s  theory  by  including  the  propeller  rotaiional 
speed  in  the  correlation  area,  a  term  defined  in  [1].  His  results  show  only 
small  humps  at  blade  rates  and  were  again  unable  to  explain  the  near 
blade-rate  humps  demonstrated  in  experiments. 

In  the  present  report,  the  correlation  analysis  follows  Sevik’s  approach 
[1],  but  removes  the  simplification  of  the  velocity  correlation  tensor  used 
in  [1]  and  [8].  The  blade  rotational  effect  is  correctly  preserved  in  the  anal¬ 
ysis  and  the  results  resemble  the  experimental  observation.  Numerically, 
for  a  given  turbulence  correlation  function,  the  statistical  properties  of  the 
unsteady  angle  of  attack  encountered  by  the  moving  blade  are  evaluated 
at  each  time  step.  Then  the  force  spectrum  is  obtained  from  a  Fourier 
transformation  of  the  computed  time  history  of  those  unsteady  angles  of 
attack.  The  next  two  sections  outline  the  fundamental  approaches  and 
the  corresponding  computational  procedures.  Numerical  examples  are  in¬ 
cluded  for  the  thrust  spectrum  of  isotropic  turbulence  on  a  rotor  using 
both  correlation  and  spectrum  approaches.  These  results  are  compared 
with  existing  experiments.  A  flow  chart  of  the  present  numerical  tur¬ 
bulence  ingestion  (TI)  model  is  given  in  Figure  1;  the  blade  surface  is 
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divided  into  chordwise  strips  along  the  span.  Each  strip  is  considered  as 
a  two-dimensional  (2-D)  section  with  unsteady  lift  concentrated  at  the 
quarter-chord  of  that  section.  As  is  indicated  in  Figure  1,  the  two  criti¬ 
cal  components  of  the  model  are  the  inflow  turbulence  characteristics  and 
the  transfer  function  between  the  turbulence  intensity  spectrum  and  the 
unsteady  forces.  In  the  given  examples,  the  Sears  function  is  used  as  the 
fluid  dynamic  transfer  function. 


Theoretical  Approach 

The  analysis  of  turbulence- generated  propeller  broadband  force  can  be 
approached  with  two  different  techniques.  One  is  the  spectrum  approach 
and  the  other  is  the  correlation  approach.  These  two  approaches  should 
lead  to  the  same  results  if  the  specified  spectrum  form  and  correlation 
function  describe  the  same  turbulence  characteristics.  The  authors  prefer 
the  correlation  method  because  it  is  the  more  direct  approach  and  is  easier 
to  apply  numerically.  However,  for  completeness  of  the  study,  the  spec¬ 
trum  approach  is  included  in  this  report.  Results  from  both  approaches 
will  be  compared  in  the  discussions. 

Spectrum  Approach 

This  section  follows  the  procedure  described  in  Blake  [7],  which  assumes 
that  the  turbulence  is  convected  frozen  (Taylor’s  hypothesis)  and  the  tur¬ 
bulent  characteristics  are  represented  as  a  product  of  functions  in  each 
direction  (separation  of  variables).  With  the  turbulent  inflow  component 
normal  to  the  blade,  U2  in  Figure  2,  expressed  in  terms  of  its  Fourier 
transform: 


(1) 


where  ^  denotes  the  vector  in  terms  of  the  coordinates  in  the  propeller 

plane  (a:i,a:3)  and  its  normal  (2:2),  the  lift  per  unit  radius  on  the  s-th 

blade  written  in  terms  of  the  Fourier  components  of  the  incident  upwash, 
— ♦ 

U2{k,uj),  is 

dL.(u)/dr  =  p7rC«2(t,a.0!7(r)L2p(/:iC/2)e'(*»’-+‘-'=''*-"''  (2) 

where  L2D  Is  the  two-dimensional  lift  due  to  a  unit  disturbance,  C  is 
the  blade  chord  length,  b  is  the  blade  spacing  at  the  radius  r,  A-i2  equals 
A:i  cos4>  -f  A:2sin<;S',  (}>  is  the  blade  pitch  angle,  and  U(r)  is  the  resultant 
velocity  into  the  blade  at  radius  r.  The  net  force  on  an  N-bladed  propeller 
in  the  i-th  direction  due  to  the  turbulence  contribution  is  found 

by  summing  over  all  the  blades  and  integrating  over  the  radius  from  the 
propeller  hub  (Rh)  to  tip  (Rt) 

F.ik.u)  =  /""p7rCu2(t,«)i7(r)ni(.ii)i2o(iiC/2)  (3) 

s  =  0 

where  n,(0)  is  the  direction  cosine  in  the  i-th  direction  of  the  blade.  The 
summation  over  the  N  blades  in  equation  (  3)  is  called  a  filtering  function, 
denoted  and  is 

.4,(A:i26)  - 

5=0 

sin(Aaq2V2)et»(^-bb26/2] 

sin(A’i26/2) 

Substituting  this  expression  back  into  equation  (  3),  the  net  force  becomes 
F,{lu;)  =  /"  p-Ci-2(e-')C(r)7,,(«i)L20(i-,C/2).4,(i-,2i)e-'^”— "rfr  (5) 

J  tin 


(4) 


and  the  frequency  spectrum  of  the  fluctuating  force,  becomes 

=  i  i  i  C  C  I  1^1  Mh.b) 

U{ri)U{r2)^22{k,uj)ni{-fi)ni{^2)^'‘'''^''^~'''^  dr2dridk  (6) 

where  <E>22(^?  w)  is  the  wave  number  spectrum  of  the  turbulent  field.  Equa¬ 
tion  (  6)  is  the  fundamental  form  of  the  spectrum  approach.  It  is  clearly 
seen  that  the  major  task  of  the  analysis  is  in  the  specification  of  the  wave 
number  spectrum  of  the  turbulent  field  ^-^d  the  evaluation  of 

the  multiple  integrals.  In  general,  the  spectrum  form  $22(^»^)  associated 
with  a  moving  body  is  not  known,  and  the  prediction  of  'I'„(ixj)  with  triple 
integration  over  k  is  impractiable.  Neverthcness,  the  analytical  results 
associated  with  some  simple  forms  of  $22(^>‘^)  could  provide  a  good  ap¬ 
proximation  to  reality  and  lead  to  a  better  understanding  of  the  overall 
features  of 

By  applying  Taylor’s  hypothesis  that  the  turbulence  convects  over  the 
blade  with  the  resultant  velocity  and  without  decay,  and  by  assuming 
that  the  wave  number  spectrum  can  be  separated  into  directions  along  the 
resultant  velocity  (^i),  the  radial  direction  (A:3),  and  the  normal  direction 
(A:2),  the  turbulence  spectrum  becomes 


<I>22(^,^’)  =  ul4>i{ki)<f)2ik2)(f>z{k:i)8[u  -  (7) 


Furthermore,  we  assume  that  takes  a  form  of 


Mki) 


1  A. 

-l  +  {Ark,y 


(8) 
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where  2At  represents  the  integral  length  scale  of  U2  in  the  i-th  wave  vector 
direction.  When  the  radial  integral  length  scale  is  small  with  respect  to 
the  blade  span,  i.e., <^3(^:3)  ~  A3/X  and  the  radial  direction  wave  number 
is  small,  i.e.,  k:iR  <C  1  ,  the  force  spectrum  ,  equation  (  6)  becomes 

=  Jk,  k  C  I  L,D(kiC/2)  PI  A.(k:2h)  p 

U‘^{r)u2(f>i(ki)(l)2ik2)n‘f{(f))S[ijj  —  U{r)ki]drdkidk2 
""  I  A,[(/:2sin<?i>  +  ^cos<^)6]  ^ 

U‘^{r)u\4'\{u /U)(t>2{k2)'n]{(p)drdk2  (9) 

For  simplification,  the  average  spanwise  lift  amplitude  as  a  function  of 
wave  number  has  been  used  in  the  last  formulation.  Then  the  turbu¬ 
lence  spectrum  can  be  moved  outside  of  the  radial  integral.  The  discrete, 
numerical  expression  of  equation  (  9)  becomes 

=  L  <f>2{k2)  E  A3  I  L2d{^)  Pi  As[{k2  sin  (/>  -I-  ^  cos  <f))b]  p 

6k2  SR  ^ 

U'^ul<f)i(ujlU)n-{<f))  (10) 

Equation  (  10)  entails  summation  over  the  blade  elements  and  the  wave 
number  and  it  is  numerically  more  managable  than  :quation  (  6).  The 
application  of  this  model  will  be  presented  in  a  later  section  and  the  results 
will  be  discussed.  For  the  case  of  small  advanced  coefficient  and  small  hub, 
Rt  Rhi  equation  (  9)  can  be  further  reduced  to 


I  £2d(^)  P  (^I^i/Rt) 

"<('>')^'’^‘(^)£1  I  ^s[(*^2sin7  +  ^cos7)6]  p 

4>2{k2)dk2  (11) 


where  Ut  is  the  resultant  tip  speed  sjV^  +  {RtO.)^  (V  is  axial  inflow  ve¬ 
locity).  Blake[7]  obtained  the  asymptotic  solutions  of  the  above  equation. 
Evaluation  of  the  integral  is  determined  by  the  relative  width  of  the  turbu¬ 
lence  spectrum,  <^2(^2),  and  the  filter  function,  |  A5[(  A;2  sin  7  +  ^  cos  7)6]  p. 
When  the  band  of  <^2(^2)  is  greater  than  the  range  of  the  filter  function, 
i.e.,  the  axial  integral  length  parameter,  2A2,  is  smaller  than  the  blade 
spacing  projected  in  the  axial  direction  (small  correlation  length),  the  last 
equation  can  be  simplified  to 


— rinV2N2D  4/  C  N2/7r.2^2  2A3 

^  1  +  {Aiuj/asRry 

1 

1  -}-  t:uC /^sRt 

(12) 


On  the  other  hand,  when  the  axial  integral  length  parameter,  2A2,  is 
greater  than  the  blade  spacing  projected  in  the  axial  direction  (large  cor¬ 
relation  length),  equation  (11)  becomes 


8A’'x2  1  2>.2  n  4  C  .2/^  n2  ^2  2-^3  Ai/ttRt 

3ns  ^  ^Rt^  Rt1  +  (AimA7i?r)2 

I  ^[(^2  sin  7 -I- ^  cos  7)6]  p 
1  ttuiNC/Rt 
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where  m  is  the  number  of  the  blade  rate  frequency;  otherwise, 
equals  zero. 

In  examining  equations  (  12)  and  (  13),  it  is  seen  that  equation  (  12)  is 
a  monotonically  decreasing  function  of  uj  and  it  resembles  the  rotational 
uncorrelated  results  of  Sevik  [1].  Equation  (  13)  contains  only  the  blade 
rate  humps  with  non-decayed  amplitude  above  the  rotational  uncorre¬ 
lated  results.  Neither  of  these  approximations  can  be  applied  to  predict 
the  complete  measured  phenomena.  In  order  to  obtain  a  meaningful  pre¬ 
diction,  one  has  to  go  back  to  at  least  equation  (  11).  Numerical  results 
show  that  even  (11)  does  not  give  satisfying  predictions;  only  the  results 
of  equation  (  10)  are  in  reasonable  agreement  with  the  measurements. 
Further  discuissions  of  the  spectrum  approach  will  be  given  in  the  results 
section. 


Correlation  Approach 

The  derivations  of  the  frequency  spectrum  of  propeller  thrust  using  the 
correlation  approach  are  briefly  discussed  here.  These  derivations  follow 
the  theory  of  Sevik  [1]  except  that  the  velocity  correlation  has  been  mod- 
ifled  to  incorporate  propeller  rotational  effects.  In  order  to  calculate  the 
characteristics  of  the  fluctuating  force  over  the  blade  surfaces,  the  blade 
is  divided  into  a  number  of  surface  elements.  The  time-dependent  forces 
acting  on  the  various  surface  elements  are  related  by  virtue  of  spatial  and 
temporal  correlation  of  the  velocity  fluctuations  in  the  approach  stream 
as  well  as  by  virtue  of  the  induced  effects  that  take  place  between  adja¬ 
cent  elements.  In  the  following  tensor  equations,  subscripts  are  used  to 
denote  the  direction  along  the  coordinate  axes,  while  superscripts  are  used 
to  denote  the  blade  element  involved.  For  example,  i^j^(f,T)  denotes  the 
hydrodynamic  force  acting  on  the  a-th  element  in  the  direction  i  at  the 
instant  of  time  t  caused  by  a  velocity  fluctuation  of  unit  magnitude  in  the 
direction  j  on  the  ^-th.  element  at  the  instant  r.  With  this  convention,  the 
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hydrodynamic  force  acting  on  the  o-th  element  at  time  t  in  the  direction 
i,  due  to  the  influence  of  velocity  fluctuation  at  all  elements,  /?,  and 

in  all  the  directions,  over  all  the  time,  r,  is  expressed  as 

=  f^Ft/{t,r)u^{T)dT  (14) 

where 


=  1>2,3 
a,l3  =  l,2,...n 


Since  lf{t)  is  a  random  function  of  time  in  a  turbulent  flow,  a  statistical 
approach  must  be  employed.  If  the  inflow  turbulence  is  time-invariant,  the 
correlation  function  of  unsteady  forces,  becomes 

<(t)  =  Em)t^{t  +  r)] 

=  f  r  r  r  > 

ul^{t  T  —  T2)dTidT2dt 

=  +  Ti  -  T2)dTidT2  (15) 

where 

-  ^2)  =  ^  u]{t  -  ri)u^(f  +  T-  T2)dt 

The  corresponding  frequency  spectrum  of  the  correlated  force  fluctua¬ 
tions  on  the  o-th  element  in  the  direction  i  due  to  the  /3-th  element  in 
the  direction  j  can  be  found  by  taking  the  Fourier  transform  of  the  above 
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correlation  tensor: 


<(-)  = 

=  ^  Pik('r\)driFfi(Ti)dTiRf„(T +  Ti  -  T2)dT 

=  \Hri;'{u^)nHfi(u)]Gt(u)  (le) 

where 

=  ffT{r)e-^'dT 

and 

H°‘if{u})  is  the  hydrodynamic  frequency  response  function  and  Otmiu))  is 
the  Fourier  transform  of  the  velocity  correlation  function.  Equation  (  16) 
is  the  formula  for  predictions  of  unsteady  forces  due  to  inflow  turbulence. 
When  we  apply  this  equation  to  the  rotors,  we  have  to  specify 
and  according  to  the  enviroments  in  which  the  rotors  are  operat¬ 

ing.  A  computer  code  has  been  written  to  compute  the  G(u;)  and  for 
specified  turbulence  characteristics.  To  demonstrate  the  computation  pro¬ 
cedure,  the  analysis  of  an  isotropic  turbulence  is  given  below  for  general 
discussion. 


Isotropic  Turbulence 

In  the  present  method,  we  shall  first  derive  the  velocity  correlation 
function  between  points  o  and  A  detailed  derivation  for  isotropic  tur¬ 
bulence  is  given  by  Hinze  [9]  and  is  briefly  described  below.  Longitudinal 
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correlation  refers  to  the  coordinate  system  in  Figure  3.  The  origin  of  the 
coordinate  system  is  at  o  and  r  is  the  distance  between  points  o  and 
with  components  along  the  axis  rr,-.  In  the  plane  of  the  o /3-line  and 
the  a:i-axis,  the  velocities  at  a  and  /3  can  be  resolved  into  the  components 
ui*,  directed  between  these  two  points,  and  the  components  which  are 
perpendicular  to  the  line-o^.  The  velocity  components  along  the  s’! -axis 
at  these  two  points  are: 


=  u®*  cose,:?!  —  U2»  sin 


Let  us  define 

^2*^2*  = 

where  u  is  the  root-mean-square  of  turbulence  and  r  is  the  distance  along 
these  two  points.  Because  of  invariance  conditions,  we  have 

When  we  multiply  u®  by  itf  and  take  the  mean  value  of  the  product,  we 
get 


12 


i2ii  =  UiUi 


=  ^  [ - ;:2 - 


66  +  p(^)] 


(18) 


Next,  let  us  consider  From  Figure  4,  we  have 


u 


a 

1 


U 


.6 


u 


if 


(19) 


where  F{u^^)  is  the  term  associated  with  the  velocity  component  uf,.  This 
term  is  not  given  in  detail  because  the  invariance  and  isotropy  conditions 
provide  the  zero  mean  for  all  terms  containing  the  uf*  component.  The 
invariance  and  isotropy  conditions  yield 


= 

=  U2»Ui^ 


_  ,,o  ,,,0 

=  0 


The  correlation  of  and  U2  is  determined  as 


R 
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W?W2 
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(20) 


=  u 


2/(^)  -  9{r) 


66 


In  general,  combining  equations  (  18)  and  (  20),  we  can  write  the  ve¬ 
locity  correlation  as 


Rii  =  (21) 

where  6,j  is  the  Kronecker  delta.  The  relationship  between  /(r)  and  g{r) 
can  be  derived  from  the  continuity  condition  of  fluid  incompressibility. 


dRij 


=  0 


Substituting  Rij  of  equation  (21)  into  equation  (  22)  yields 


/(^)  + 


rdfjr) 
2  dr 


9{r) 


(23) 


From  this  relation  between  f{r)  and  g(r),  the  velocity  correlation  function 
can  be  expressed  in  terms  of  one  scalar  function,  either  f(r)  or  g(r),  which 
gives 


2  dr 


(24) 


The  longitudinal  correlation  function,  /(r),  can  be  approximated  by  an 
exponential  function 
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(25) 


where  A  is  the  integral  scale  of  the  turbulence.  This  is  equivalent  to 
assuming  that  the  turbulence  spectrum  is  in  the  form  as  it  is  used  in  the 
spectrum  approach.  As  has  been  verified  by  screen-generated  turbulence, 
this  is  a  good  approximation  for  homogeneous  and  isotropic  flow  field. 

For  blades  rotating  in  a  turbulent  field  with  rotational  speed  r(r)  is 


r(r)  =  )J{VTy  +  rl  +  -  2rar0cos{9a  -  00  +  fir)  (26) 


Equation  (  24)  can  be  rewritten  as: 


(27) 


or  in  the  algebraic  form 


Ru 
R22 
^12  = 


r  _nil 
— )  e  A 


Ofl  O  / 

"  -  5a""'  ^ 


_Ill2 


^2rA 


(28) 


By  taking  numerical  Fourier  transformations  of  the  time  domain  cor¬ 
relation  function  (equation  (  Tl')'),  oi^e  ih>tams  Sh\eVrequencs  'specvra'm 
the  velocity  correlation  function^  Gij{u),  as 
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GiiH  =  i  /_"  %(r)e-'“Vr  (29) 

where  the  omitted  superscripts  axe  the  indices  of  two  points  for  conve¬ 
nience. 

In  examining  equation  (  26),  it  is  clear  that  Rij  is  periodic  with  respect 
to  blade  rotation  due  to  the  presence  of  the  cosine  term  in  the  correlation 
function.  This  periodic  term  is  most  dominant  when  V/Q  is  small.  As  r  or 
V  increases,  the  importance  of  the  cosine  term  decreases  and  Rij  becomes 
a  monotonically  decreasing  function  of  r.  The  spectrum  represented  by 
equation  (  29)  should  thus  contain  both  a  broadband  and  a  blade  rate 
hump  spectrum  in  general.  It  resembles  the  experimental  measurements 
and  will  be  discussed  in  detail  in  the  results  section. 

Sevik  [1]  and  Martinez  [8]  studied  special  cases  of  the  above  equation. 
Sevik  [1]  considered  the  case  where  Q.R  is  much  smaller  than  V  and  he 
approximated  equation  (  26)  by 


r{T  =  0)  =  tJrl  +  Tj-  2Tar0  COs(Sa  -  0/)} 


(30) 


The  corresponding  velocity  correlation  then  becomes 


r(T=0) 

RijC:^e  A  R^i{VT) 


(31) 


The  components  of  the  velocity  correlation  tensor  are 


Rii  =  u  e 


2.-(¥+-^) 


(32) 
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i?22  W  (1 


o  ,  Vt  .  /Vt  ,  r(T  =  0)^ 
2/1 - V~(  — +  A') 


The  effect  of  the  angular  location  of  the  blade,  fir,  is  totally  absent  in 
Sevik’s  correlation  model;  is  a  time  independent  term.  The  so¬ 

lution  for  Gij  contains  no  blade  humps  and  may  be  expressed  analytically; 
that  is, 


9 

— - — - e" 

X  (^)2+u;2 

2  ^ 

-u‘^ _ A _ (I  _  - 


1  /£\2  _  ,  ,2 
2(^)2 +u>2' 


Hoping  to  recover  the  blade  rate  humps  without  losing  the  simplicity 
of  the  analysis,  Martinez  [8]  modified  Sevik’s  analysis  by  including  the 
propeller  rotational  speed  through  the  approximation 


-  V^Q  +  ’^0-  cos{9a  -00-]- 


(35) 


Then  the  velocity  correlation  function  becomes 


2  -f— +^1 

—  V  A  ^  A  / 


Rii  =  u  e 


(36) 


Jin  =  nil-  —)e 


(37) 
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t 

which  still  preserves  the  exponential  form  of  the  correlation  function  used 
by  Sevik  [l].  The  analytical  form  of  the  corresponding  spectrum  is  given 
by  Martinez  in  reference  [S|.  The  solution  contains  Bessel  functions.  The 
modification  of  Martinez  made  the  solution  more  complex  and  still  did 
not  recover  the  blade  rate  humps  as  he  expected.  As  is  true  for  the 
spectrum  approach,  simplification  of  the  correlation  function  implies  a 
filtering  process  in  the  physics.  For  practical  applications,  the  search  for 
an  analytical  solution  could  have  prevented  the  advancement  of  unsteady 
force  prediction  methods  in  the  past.  The  present  analysis  adapts  both 
theoretical  and  numerical  techniques.  For  individual  cases,  the  theoretical 
analysis  will  be  carried  to  its  full  extent,  in  this  case  to  equation  (  29), 
and  then  numerical  analysis  will  follow. 

Fluid  Dynamic  Response  Function  -  H(uj) 

The  next  step  of  the  computation  is  the  establishment  of  a  transfer 
function,  defined  in  equation  (  16).  The  simplest  computation 

is  based  on  two-dimensional,  incompressible  and  inviscid  aerodynamics 
theory.  Sears  [10]  derived  the  response  function  for  a  zero  skew,  discrete 
sinusoidal  gust  convected  in  the  plane  of  a  foil  at  zero  angle  of  attack.  The 
resulting  aerodynamic  unsteady  force  is  represented  as 

C 

H(u;)  =  2'KpU—K{u)  (38) 

where  A'(u;)  is  the  Sears’  function  and  is  in  terms  of  cylindrical  functions 
and 
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(39) 


K{u) 


H?\k) 


where  k  is  the  normalized  reduced  frequency  {u)Cf2U). 

Equation  (  38)  was  derived  for  an  isolated  foil  operating  with  sinusoidal 
velocity  disturbances  normal  to  the  foil  surface.  When  applying  the  equa¬ 
tions  to  propellers  or  small  aspect  ratio  foils,  the  results  of  using  a  two- 
dimensional  approximation  may  not  valid.  However,  the  two-dimensional 
theory  will  give  satisfactory  predictions  providing  the  value  of  the  local 
lift  slope  is  appropriately  corrected.  This  means  that  the  coefficient  277  of 
equation  (  38)  be  replaced  by  the  steady  state  local  lift  slope  as  described 
by  Sevik  [1]. 

Applying  the  above  Sears’  function  to  a  two-dimensional  blade,  the 
angle  between  the  inflow  and  the  direction  of  advance  also  has  to  be  con¬ 
sidered.  This  gives 


Hp{  uj)  =  H{uj)cos(f) 


(40) 


where  Hp  denotes  the  two-dimensional  response  function  for  a  propeller 
in  its  direction  of  advance  and  <f>  is  the  inflow  angle  with  respect  to  that 
direction.  When  operating  near  design  conditions,  (j)  may  be  approximated 
by  the  pitch  angles  of  blade  sections  for  most  marine  propellers. 


Numerical  Procedures 

With  a  specifled  correlation  function,  equation  (28),  and  a  chosen  trans¬ 
fer  function,  equation  (  40),  we  may  proceed  to  perform  the  numerical 
unsteady  force  computations  as  given  in  equation  (  16).  Since  propeller 
unsteady  forces  are  generated  from  the  unsteady  velocity  relative  to  the 
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blades,  the  velocity  spectrum  componeiits  G{co)  in  equation  (  16)  should 
be  expressed  in  a  rotational  cylindrical  coordinate  system  which  is  fixed 
on  the  blade.  Let  the  mean  inflow  to  the  propeller  be  denoted  by  i/“. 
Then  the  unsteady  angle  of  attack  becomes 


with 


(41) 


^  A’  —  <^08  6q  -f  U0  sin  <t)a  (42) 


and 


=  {u^{0)  cos  (f)Q{0)  +  U0  s'm(f>a{0))  +  {u^  cos  61^(7 )  + Uq  sm(i)0{r)) 
—  Ujuf  cos  ^q(O)  cos  60(7)  +  cos  <f)Q{0)  sin  <ji>/?(r)  + 

sin  <?5»q(0)  cos  <?i/?(T)  +  sin  (25*0(0)  sin  <f)0{7)  (43) 


where  (t)a  is  the  propeller  inflow  angle  at  point  q  and  is  defined  in  equation 
(  39). 

The  velocity  correlation  of  equation  (  16),  thus,  becomes 


=  Rir  cos  <2^0 (0 )  cos  60{7)  COS  <25>a(0 )  siu  (25>;d  T )  + 

sm  daiO)  cos  60(7)  +  sin  <2)q(0)  sin  (2>;j(r )  (44) 
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The  calculation  of  is  straight  forward  as 

where  dx  =  Vr  and 

r(r)  =  \j{dxY  +  rj  +  -  2r^,r^  cos(0q(O)  -  ^/?(0)  +  Qt) 

(46) 

The  calculations  of  and  in  terms  of  r  has  to  be  further  ex¬ 
pressed  through  the  instantaneous  angular  locations  of  the  blade  elements. 
For  the  coordinate  system  shown  in  Figure  5,  the  tangential  velocity  com¬ 
ponent  is 

Uq  =  u^': os  O0{t)  —  u°  sin  6 0{t)  (47) 

and  the  velocity  correlations  and  R^^  are 

Ke  =  Ky  cos  O0{t)  -  sin  0/?(r)  (48) 

and 

=  R^y  cos  6q{0)  cos  6  i}{t)  +  R^^  sin  6a{0)  sin  6  p{t)  — 

cos  ^a(O)  sin  O0(t)  —  Ry^  sin  0q(O)  cos  O0(t)  (49) 

where 
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and 


^yy 

'■xy 


dy 


=  u^(l  -  ■~dxdy)e~^'^ , 

=  u^(l  — -^dxdz)e~^  ^ 

=  cos0q(O)  —  r/j  cos^/j(r), 


dz  =  Ta  sin^a(O)  —  r/j  sm0/j(r). 


Since  turbulence  has  been  assumed  to  be  homogeneous,  is  only  a 

function  of  r,  and  can  be  evaluated  at  0a(O)  and  0/?(r). 

Based  on  the  above  definition  of  i?iviv,  correlation  functions  are  numer¬ 
ically  evaluated  at  given  uj  values  as 

^NJv(^)  =  ^Xoo  ^^dr  (50) 

The  thrust  spectrum  tensor,  equation  (  16),  between  blade  elements  a 
and  P  is  then  calculated  as: 


(51) 


Assuming  the  unsteady  force  generated  at  each  individual  strip  is  two 
dimensional,  the  last  equation  becomes 
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=  [H‘^(u)YlH^^(u)][GlUu)] 


with  Hp  given  by  equation  (  39).  The  total  unsteady  force  spectrum  is 
obtained  from  summing  over  all  elements  of  ct  and  0.  In  this  com¬ 

putation,  the  mean  flow  angle  <()  has  to  be  estimated  before  the  unsteady 
thrust  can  be  determined.  This  angle  can  be  obtained  from  propeller  com¬ 
puter  codes.  For  near-design  operating  conditions,  the  pitch  of  the  blade 
may  be  used  for  (f)  for  most  marine  propellers,  as  mentioned  earlier. 


Results  and  Discussions 

Sevik’s  experiment  [1]  is  used  to  study  the  present  approach  and  to 
demonstrate  the  differences  between  the  present  computation  and  previ¬ 
ous  results.  The  experiment  was  conducted  in  the  1.22m  (48”)  diameter 
water  tunnel  at  Pennsylvania  State  University.  The  propeller  used  for  this 
investigation  has  ten  blades  with  a  constant  chord  length  of  2.54cm  (1”) 
and  a  radius  of  10.16cm  (4”).  The  turbulence  level  in  the  test  section 
without  the  grid  was  about  0.1%.  Two  grids  were  used  with  mesh  sizes  of 
10.16cm  (4”)  and  15.24cm  (6”),  respectively,  to  generate  the  turbulence 
level  of  about  3%  with  two  different  integral  length  scales  at  the  propeller 
plane.  The  distance  between  the  grid  and  the  propeller  plane  was  twenty 
times  the  grid  sizes.  Data  were  taken  at  a  tunnel  speed  equal  to  4.57m/sec 
(15ft/sec)  and  the  advance  ratio  was  1.22.  Power  spectral  density  of  the 
propeller  thrust  due  to  turbulence  was  measured  and  compared  with  theo- 
retial  results  of  Sevik  [1],  Figure  6.  The  measured  spectra  show  significant 
humps  at  the  first  blade  rate  frequency  while  the  humps  at  the  higher  blade 
rates  were  observable  only  for  the  15,24cm  grid  case.  The  theoretical  re¬ 
sults  of  Sevik  compared  well  with  the  measurements  with  respect  to  the 
broadband  part  of  the  spectrum  while  the  humps  were  missing.  Later, 
Blake  [7]  compared  those  data  with  his  asymptotic  solutions  of  equation 
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(  12)  and  equation  (  13),  derived  to  improve  the  prediction  on  blade  rate 
humps.  His  results  are  also  given  in  Figure  6.  The  smooth  broken  line  in 
Figure  6  shows  the  asymptotic  results  of  small  correlation  (equation  (  12)) 
and  the  dotted  line  provides  the  large  correlation  results  (equation  (  13)). 
It  is  seen  that  the  small  correlation  results  approach  Sevik’s  theory  as 
expected  since  they  are  equivalent  to  omitting  the  rotational  correlation 
from  the  analysis.  The  large  correlation  asymptotic  results  give  humps 
centered  at  multiples  of  the  blade  rate  frequency  without  the  broadband 
part  of  the  spectrum.  The  comparison  of  the  asymptotic  solutions  and 
measurement  indicate  that  the  basic  approach  is  correct  and  that  numer¬ 
ical  simplification  should  be  avoided. 

Computations  based  on  equation  (  10)  and  equation  (  11)  were  per¬ 
formed  numerically  for  the  two  grid  sizes,  and  the  results  are  given  in 
Figures  7  and  8.  The  results  of  equation  (  11),  Figure  7,  show  that  the 
center  frequencies  of  the  humps  shift  and  skew  to  higher  frequencies  when 
compared  with  Blake’s  calculations.  The  amplitudes  of  the  humps  do  not 
decay  at  higher  blade  rates  for  either  grid  size.  Also,  the  results  using 
equation  (11)  show  ditches  between  the  humps  quite  below  the  measured 
results.  Since  equation  (  11)  did  not  have  the  radial  dependence  of  the 
filtering  function,  Sears’  function  and  Taylor’s  hypothesis,  it  over  em¬ 
phasizes  the  blade  rotational  effect.  By  removing  the  small  advance  coef¬ 
ficient  assumption,  including  the  radial  dependence  on  filtering  and  Sears’ 
functions,  and  numerically  integrating  equation  (  10),  the  hump  ampli¬ 
tude  and  ditch  depth  are  reduced  and  the  hump  bandwidth  is  increased 
as  shown  in  Figure  8.  It  is  thus  demonstrated,  in  a  practical  application, 
that  the  radial  dependence  should  be  included  in  the  analysis.  Neverth- 
less,  including  this  effect  in  the  wave  number  integration  complicates  the 
computation  and  the  spectrum  approach  is  not  recommended  at  present. 

Figures  9a  and  9b  present  the  comparisons  between  the  results  of  the 
present  correlation  theory  and  the  experiment.  The  current  correlation 
theory  predicts  both  the  humps  which  were  demonstrated  in  the  experi¬ 
ment  and  the  broadband  part  of  the  spectrum.  The  broadband  part  of  the 
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spectrum  compares  well  while  the  hump  amplitudes  and  center  frequen¬ 
cies  of  higher  harmonics  do  not  agree.  The  cause  of  the  hump  discrepancy 
could  be  in  the  uncertainty  of  the  inflow  turbulence.  The  experiments 
show  the  inconsistency  in  the  results  betweenl 0.16cm  (4”)  and  15.24cm 
(6”)  grids.  The  15.24cm  (6”)  grid  results  have  humps  at  both  first  and 
second  blade  rate  frquencies  while  10.16cm  (4”)  grid  results  show  only  one 
hump  near  the  first  blade  rate.  These  differences  could  reflect  deviations 
from  the  design  homogeneous  and  isotropic  grid  turbulent  inflow.  Since 
the  characteristics  of  the  inflow  were  not  measured,  no  further  compar¬ 
isons  are  possible.  However,  from  the  general  features  of  the  comparsion 
and  the  other  computations  which  are  not  included,  the  authors  are  cer¬ 
tain  that  the  present  model  should  predict  the  spectrum  reasonably  well 
if  the  inflow  turbulence  is  given. 

The  characteristics  of  the  unsteady  force  spectrum  may  be  general¬ 
ized  from  the  study  of  isotropic  turbulence.  The  present  isotropic  results 
show  that  the  spectrum  can  be  normalized  in  terms  of  the  parameters 
V/Afi,  A/jR  and  Their  influence  on  the  resulting  unsteady  force  is 

demonstrated  in  Figures  10  and  11,  and  is  summarized  as  follows: 

1)  The  effect  of  the  inflow  turbulence  scale,  A/i2:  In  general,  the  smaller 
the  turbulence  length  scale,  the  higher  the  broadband  unsteady  spectrum. 
However,  both  amplitudes  and  bandwidths  of  the  blade  rate  humps  de¬ 
crease  with  decreasing  turbulence  scale  as  the  result  of  faster  decay  of  the 
rotational  correlations,  while  the  variation  of  the  turbulence  length  scale 
has  only  a  minimal  effect  on  the  hump  center  frequencies.  Figure  10. 

2)  The  effect  of  advance  coefficient,  F/AQ:  With  other  nondimensional 
parameters  held  constant.  Figure  11  shows  that  as  the  advance  coefficient 
decreases,  the  blade  rate  humps  become  more  pronounced.  The  higher 
blade  rate  humps  become  significant  at  low  advance  coefficients  as  shown 
by  Blake[7]  and  discussed  previously.  In  the  normal  operating  condition 
(F/AQ  greater  than  0.5),  only  the  first  blade  rate  hump  is  expected. 

3)  Skew  of  the  hump  center  frequencies:  The  center  frequencies  of  the 
humps  shown  in  Figure  12  are  not  at  blade  rate  frequencies  and  are  skewed 
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to  the  right.  The  skew  is  caused  by  the  geometrical  pitch  of  the  blade. 
The  larger  the  pitch,  the  greater  the  skew  will  be.  The  first  order  estimate 
of  this  skewness  may  be  approximated  from  equations  (  11)  or  (  48). 

4)  Inverse  cubic  dependency  on  the  high  frequencies:  To  demostrate 
the  spectrum  behavior  at  high  frequencies,  the  theoretical  results  of  the 
grid  turbulence  for  10.16  cm  and  15.24  cm  screens  are  re-plotted  in  nor¬ 
malized  form  in  Figure  13.  The  frequency  is  normalized  by  the  blade  rate 
frequency  and  the  thrust  spectrum  is  normalized  by  the  third  power  of 
the  velocity.  At  high  frequencies,  past  the  second  hump,  the  slope  of  these 
normalized  curves  approaches  minus  three  as  expected  from  equation  (  12) 
or  equations  (34)  and  (38)  for  a  given  advance  coefficient  and  turbulence 
level.  The  two  grids  used  in  the  experiment  did  not  cover  enough  of  a 
range  in  turbulence  lengths  to  make  as  substantial  alterations  near  the 
hump  as  the  theoretical  calculation  provided  in  Figure  10. 

Because  ships  normally  operate  at  a  nearly  constant  advance  coefficient, 
it  is  of  interest  to  present  the  thrust  spectrum  in  dimensional  form  as  a 
function  of  ship  speed  or  p^'opeller  rpm  with  a  given  advance  coefficient. 
The  behavior  of  the  thrust  spectra  at  different  speeds  for  the  10.16  cm  grid 
turbulence  is  shown  in  Figure  14.  It  is  seen  that  the  hump  bandwidth  at 
the  first  blade  rate  frequency  increases  as  ship  speed  increases.  In  practice, 
one  should  anticipate  a  sharp  hump  during  lower  speed  operation. 


Conclusions 


Theoretical  predictions  of  broadband  thrust  forces  in  a  turbulent  flow 
are  presented.  The  calculations  are  carried  out  using  a  correlation  method 
and  a  spectrum  method.  The  results  of  the  correlation  theory  provide  a 
better  correlation  with  experiment  than  the  spectrum  theory,  computed 
from  an  assumption  of  homogeneous  and  isotropic  turbulence.  Turbulence 
velocity  correlation  can  be  considered  as  the  combination  of  transverse 
and  rotational  correlation.  The  effect  of  transverse  correlation  produces 


26 


an  unsteady  thrust  proportional  to  the  third  power  of  ship  speed.  The 
rotational  correlation  is  inversely  proportional  to  a  rotational  parameter 
given  by  V f  KQ! .  The  rotational  correlation  distorts  the  power  rela¬ 
tion  and  causes  the  unsteady  thrust  spectra  to  show  humps  near  the  first 
and  second  blade  rate  frequencies;  these  humps  skew  to  the  higher  fre¬ 
quency  side  of  the  blade  rate  frequencies.  Since  the  hump  phenomena 
were  not  explained  previously,  only  limited,  inconclusive  experiments  are 
available  to  verify  the  theory.  Future  experiments  should  include  simulta¬ 
neous  measurement  of  unsteady  force  and  turbulence  in  order  to  validate 
the  theory. 

The  spectrum  method  is  not  as  easy  to  apply  as  the  correlation  method 
because  it  involves  higher  order  integration.  When  the  integration  or¬ 
der  was  reduced  through  simplification,  the  results  were  also  degraded  as 
shown  in  Figures  6  through  8.  The  correlation  method  gives  good  results, 
Figure  9,  and  can  be  adapted  to  more  complex  cases  with  little  numerically 
difficulty. 

In  addition  to  improving  the  inflow  turbulence  assumptions,  the  analyt¬ 
ical  methods  can  be  further  improved  by  including  the  effect  of  finite  span, 
camber  and  angle  of  attack.  A  parametric  investigation  of  the  turbulence 
inflow  and  propeller  unsteady  response  should  be  done  under  laboratory 
conditions  before  the  theory  is  applied  to  a  complex  propeller /rotor  ge¬ 
ometry. 
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Propeller  Blade  Geometry 


Propeller  Inflow  Velocity 
& 

Turbulence  Characteristics 
(isotropic,  homogeneous) 


Fluid  Dynamic  Transfer  Function 


Fig.1.  Flow  chart  of  excitation  force  model. 
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Fig. 6a  10.16cm  Grid 
Fig. 6  Comparison  of  Existing  Methods 
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Fig. 6b  15.24cm  Grid 
Fig. 6  (Continued) 


Fig. 7a  10.16cm  Grid 

Fig. 7  Comparison  of  Simplified  Spectrum  Approaches 


Present  -  Spectrum  Eq.ll 
Blake  -  small  correlation 
Blake  -  large  correlation 
Experiment  [1] 
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Fig.  7b  15.24cm  Grid 
Fig. 7  (Continued) 


Present  -  Spectrum  Eq  10 
Present  -  Spectrum  Eq  ll 
Experiment  [1] 
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Fig. 8a  10.16cm  Grid 

Fig. 8  Comparison  of  Unified  Spectrum  Methods  with  Experiment 
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Fig. 8b  15.24cm  Grid 
Fig. 8  (Continued) 
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Fig. 9a  10.16cm  Grid 

Fig. 9  Comparison  of  Correlation  Methods  With  Experiment 
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g.9b  15.24cm  Grid 
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Frequency,  Hz 

Fig.  10  Effect  of  Turbulence  Length  Scale  (A/R)  on  Unsteady  Thrust 

(Turbulence  Level  =  0.03;  V/RQ=:0.39) 
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Frequency,  Hz 

Fig.  11  Effect  of  Advance  Coefficient  (V/AH)  on  Unsteady  Thrust 

(Turbulence  Level  =  0.03;  A/R=0.28) 
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Fig.  12  Normalized  Effect  of  Advance  Coefficient  (V/Afi)  on  Unsteady  Thrust 

(Turbulence  Level  =  0.03;  A/R=0.28) 


Normalized  Frequency,  w/fi' 

[3  Normalized  Unsteady  Thrust 
(Turbulence  Level  =  0.03;  J=1.2^ 


Fig.  14  Unsteady  Thrust  at  Constant  Advance  Coefficient 
(Turbulence  Level  =  0.03;  10.16cm  Grid,  J=1.22) 
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